/*******************************************************************************

This code file produces Figures 4, A14, and A15.

*******************************************************************************/

*** Manage settings

	run "$dir/code/modules/settings.do"
	
*** Load data

	use "$data/clean/cleaned_data.dta", clear

********************************************************************************
* Store estimates of baseline models for subsequent structural analyses
********************************************************************************

	* Estimate logit model, w/ FE
	
	logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.nta i.yearpermit, cl(nta) difficult
	mat b = e(b)
	local coef_dtaxrate = b[1,1]
	
	* Simulate 0.1 p.p. tax change for on-site inclusionary units
	predict xb, xb
	gen pr = exp(xb)/(1+exp(xb))
	gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
	
	tempfile tmp_fe
	save `tmp_fe', replace
	
	* Estimate logit model, w/o FE
	
	use "$data/clean/cleaned_data.dta", clear
	
	logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.borough i.yearpermit, cl(nta)
	mat b = e(b)
	local coef_dtaxrate = b[1,1]
	
	* Simulate 0.1 p.p. tax change for on-site inclusionary units
	predict xb, xb
	gen pr = exp(xb)/(1+exp(xb))
	gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
	
	tempfile tmp_nofe
	save `tmp_nofe', replace
		
********************************************************************************
* Structural analyses, by neighborhood
********************************************************************************	
	
*** Elasticity

	** With NTA FE

		use `tmp_fe', clear
		
		* Collapse to neighborhoods
		collapse (mean) pr pr_plus [aw=unitsres], by(nta_code)
		
		* Calculate elasticity
		gen elasticity = 1000*(pr_plus - pr)
		
		* Statistics on neighborhood-level elasticity distribution	
		summ elasticity
		kdens elasticity
			
		* Map elasticity by neighborhood
		
		tempfile elasticity_fe
		save `elasticity_fe', replace
		
		use "$data/GIS_boundaries/neighborhood_boundaries/df.dta", clear
		merge 1:1 nta_code using `elasticity_fe'
		
		format elasticity %6.2f

		spmap elasticity using "$data/GIS_boundaries/neighborhood_boundaries/coords.dta", id(_ID) fcolor(Greys2) cln(5) ndf(white) legend(position(4)) osize(vthin vthin vthin vthin vthin) ndsize(vthin)
		
		graph export "$figs/elasticity_nta_fe.pdf", replace
		graph export "$figs_overleaf/elasticity_nta_fe.pdf", replace
		
	** With Borough FE

		use `tmp_nofe', clear
		
		* Collapse to neighborhoods
		collapse (mean) pr pr_plus [aw=unitsres], by(nta_code)
		
		* Calculate elasticity
		gen elasticity = 1000*(pr_plus - pr)
		
		* Statistics on neighborhood-level elasticity distribution	
		summ elasticity
		kdens elasticity
			
		* Map elasticity by neighborhood
		
		tempfile elasticity_nofe
		save `elasticity_nofe', replace
		
		use "$data/GIS_boundaries/neighborhood_boundaries/df.dta", clear
		merge 1:1 nta_code using `elasticity_nofe'
		
		format elasticity %6.2f

		spmap elasticity using "$data/GIS_boundaries/neighborhood_boundaries/coords.dta", id(_ID) fcolor(Greys2) cln(5) ndf(white) legend(position(4)) osize(vthin vthin vthin vthin vthin) ndsize(vthin)
		
		graph export "$figs/elasticity_nta_nofe.pdf", replace
		graph export "$figs_overleaf/elasticity_nta_nofe.pdf", replace

*** Marginal fiscal cost

	** With NTA FE
			
		use `tmp_fe', clear
		
		gen assessprop = pr * ((assesstot - assessland)/sh_taxable) * dtaxrate_onsite
		gen assessprop_plus = pr_plus * ((assesstot - assessland)/sh_taxable) * (dtaxrate_onsite + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = `sh_inclusionary' * pr * unitsres 
		gen inclusionaryunits_plus = `sh_inclusionary' * pr_plus * unitsres
		
		collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus (mean) xcoord_latlong ycoord_latlong, by(nta_code)
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)

		* Statistics on neighborhood-level marginal fiscal costs
		summ mcostperunit
		kdens mcostperunit
		
		* Map marginal fiscal cost by neighborhood
		
		tempfile mcost_nta_fe
		save `mcost_nta_fe', replace
		
		use "$data/GIS_boundaries/neighborhood_boundaries/df.dta", clear
		merge 1:1 nta_code using `mcost_nta_fe'
		
		replace mcost = mcost/1000
		format mcost %6.0f

		spmap mcostperunit using "$data/GIS_boundaries/neighborhood_boundaries/coords.dta", id(_ID) fcolor(Greys2) cln(5) ndf(white) legend(position(4)) osize(vthin vthin vthin vthin vthin) ndsize(vthin)
		
		graph export "$figs/mcost_nta_fe.pdf", replace
		graph export "$figs_overleaf/mcost_nta_fe.pdf", replace
		
	** With Borough FE 

		use `tmp_nofe', clear
		
		gen assessprop = pr * (assesstot - assessland) * dtaxrate_onsite/sh_taxable
		gen assessprop_plus = pr_plus * (assesstot - assessland) * (dtaxrate_onsite/sh_taxable + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = `sh_inclusionary' * pr * unitsres 
		gen inclusionaryunits_plus = `sh_inclusionary' * pr_plus * unitsres
		
		collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus (mean) xcoord_latlong ycoord_latlong, by(nta_code nta)
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)

		tempfile mcost_nta_nofe
		save `mcost_nta_nofe', replace
		
		use "$data/GIS_boundaries/neighborhood_boundaries/df.dta", clear
		merge 1:1 nta_code using `mcost_nta_nofe'
		
		replace mcost = mcost/1000
		format mcost %6.0f

		spmap mcostperunit using "$data/GIS_boundaries/neighborhood_boundaries/coords.dta", id(_ID) fcolor(Greys2) cln(5) ndf(white) legend(position(4)) osize(vthin vthin vthin vthin vthin) ndsize(vthin)
		
		graph export "$figs/mcost_nta_nofe.pdf", replace
		graph export "$figs_overleaf/mcost_nta_nofe.pdf", replace
		
		
*** Bootstrap estimation of standard errors

	** Elasticity

	tempfile tmp
	clear
	set obs 1
	gen tmp = .
	save `tmp', replace
	
	set seed 1234

	forvalues b = 1(1)1001 {
	
		di `b'
		
		quietly {
		
		use "$data/clean/cleaned_data.dta", clear
		
		bsample, cluster(nta)
		
		logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.borough i.yearpermit, cl(nta)
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		* Simulate 0.1 p.p. tax change for on-site inclusionary units
		predict xb, xb
		gen pr = exp(xb)/(1+exp(xb))
		gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
		
		collapse (mean) pr pr_plus [aw=unitsres], by(nta_code)
		
		* Calculate elasticity
		gen elasticity = 1000*(pr_plus - pr)
		
		append using `tmp'
		save `tmp', replace
		
		}
		
	}
	
	use `tmp', clear
	collapse (sd) elasticity, by(nta_code)

	merge 1:1 nta_code using "$data/GIS_boundaries/neighborhood_boundaries/df.dta", nogen
	
	format elasticity %12.2f

	spmap elasticity using "$data/GIS_boundaries/neighborhood_boundaries/coords.dta", id(_ID) fcolor(Greys2) cln(5) ndf(white) legend(position(4)) osize(vthin vthin vthin vthin vthin) ndsize(vthin)
	
	graph export "$figs/elast_nta_stderr.pdf", replace
	graph export "$figs_overleaf/elast_nta_stderr.pdf", replace
	

	** Marginal fiscal cost

	tempfile tmp
	clear
	set obs 1
	gen tmp = .
	save `tmp', replace
	
	set seed 1234

	forvalues b = 1(1)1001 {
	
		di `b'
		
		quietly {
		
		use "$data/clean/cleaned_data.dta", clear
		
		bsample, cluster(nta)
		
		logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.borough i.yearpermit, cl(nta)
		mat b = e(b)
		local coef_dtaxrate = b[1,1]
		
		* Simulate 0.1 p.p. tax change for on-site inclusionary units
		predict xb, xb
		gen pr = exp(xb)/(1+exp(xb))
		gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
		
		gen assessprop = pr * ((assesstot - assessland)/sh_taxable) * dtaxrate_onsite
		gen assessprop_plus = pr_plus * ((assesstot - assessland)/sh_taxable) * (dtaxrate_onsite + 0.001)
		
		local sh_inclusionary = 0.20
		gen inclusionaryunits = `sh_inclusionary' * pr * unitsres 
		gen inclusionaryunits_plus = `sh_inclusionary' * pr_plus * unitsres
		
		collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus (mean) xcoord_latlong ycoord_latlong, by(nta_code)
		
		gen mcostperunit = (assessprop_plus-assessprop)/(inclusionaryunits_plus-inclusionaryunits)
				
		append using `tmp'
		save `tmp', replace
		
		}
		
	}
	
	use `tmp', clear
	collapse (sd) mcostperunit, by(nta_code)

	merge 1:1 nta_code using "$data/GIS_boundaries/neighborhood_boundaries/df.dta", nogen
	
	replace mcost = mcost/1000
	format mcostperunit %6.0f

	spmap mcostperunit using "$data/GIS_boundaries/neighborhood_boundaries/coords.dta", id(_ID) fcolor(Greys2) cln(5) ndf(white) legend(position(4)) osize(vthin vthin vthin vthin vthin) ndsize(vthin)
	
	graph export "$figs/mcostperunit_nta_stderr.pdf", replace
	graph export "$figs_overleaf/mcostperunit_nta_stderr.pdf", replace

